To complement the fixed window analysis focused on emergence and elimination events, we also want to show that the dynamics of EWS reflect those of \(R_E\). I estimated \(R_E\) over the observations via particle filtering and state estimation. Here, I read in those values for each city and correlate them with top-performing EWS. Top-performing EWS are those with high AUC in the emergence and elimination scenarios. Note that \(R_E\) was calculated as:

\[ R_E = \frac{\beta_t}{\gamma} \frac{S_t}{N_t}. \]

library(tidyverse)
cities <- c("Agadez", "Maradi", "Niamey", "Zinder")  # the four city names
reffs <- tibble()  # empty tibble for storage
for(do_city in cities){
  fname <- paste0("../../results/filtered-states-", do_city, ".RDS")
  raw_filter <- readRDS(fname)
  reff_id <- which(raw_filter$state == "effective_r_seasonal")
  tmp_reff <- raw_filter$data[[reff_id]] %>%
    dplyr::select(time, med, week, date) %>%
    dplyr::mutate(city = do_city)
  
  reffs <- bind_rows(reffs, tmp_reff)
}
ggplot(reffs, aes(x = date, y = med, color = city)) +
  geom_hline(aes(yintercept = 1), linetype = 2, color = "grey35") +
  geom_line() +
  labs(x = "date", y = expression(italic(R)[E]), 
       title = "Effective reproduction number over time") +
  ggthemes::scale_color_colorblind()

Now I calculate the EWS at each observation time using the spaero::get_stats() function. I use a bandwidth of 30 weeks and set backward_only = TRUE so that only past data are used when calculating EWS.

library(spaero)
fname <- "../../data/clean-data/weekly-measles-incidence-niger-cities-clean.RDS"
measles_data <- readRDS(fname) %>%
  filter(year > 1994)  # drop first NA year, only used for modeling
all_stats <- tibble()  # empty tibble for storage
for(do_region in unique(measles_data$region)){
  
  cases <- measles_data %>%
    filter(region == do_region) %>%
    pull(cases)
  
  city_stats <- spaero::get_stats(
    x = cases,
    center_trend = "local_constant", 
    center_kernel = "uniform", 
    center_bandwidth = 30, 
    stat_trend = "local_constant", 
    stat_kernel = "uniform", 
    stat_bandwidth = 30, 
    lag = 1, 
    backward_only = TRUE
  )$stats
  
  city_stats_tb <- as_tibble(city_stats) %>%
    mutate(
      time_iter = 1:n(),
      date = unique(measles_data$date)
    ) %>%
    gather(key = ews, value = value, -time_iter, -date) %>%
    mutate(region = do_region)
  
  all_stats <- bind_rows(all_stats, city_stats_tb)
}
# best_ews <- c("autocorrelation", "autocovariance", "mean", "variance")
# all_stats <- all_stats %>%
#   dplyr::filter(ews %in% best_ews)

Now I look at scatterplots of \(R_E\) versus EWS.

reffs <- reffs %>%
  dplyr::mutate(region = paste0(city, " (City)")) %>%
  dplyr::select(-city)
ews_reffs <- left_join(all_stats, reffs)
Joining, by = c("date", "region")
ews_reffs_subcrit <- ews_reffs %>%
  filter(med < 1)
ggplot(ews_reffs_subcrit, aes(x = med, y = value, color = region)) +
  geom_point(alpha = 0.2) +
  facet_wrap(~ews, scales = "free") +
  ggthemes::scale_color_colorblind()

# for(do_region in unique(ews_reffs$region)){
#   print(ggplot(filter(ews_reffs, region == do_region & ews == "variance"), 
#        aes(x = med, y = value, color = lubridate::week(date))) +
#     geom_point() +
#     facet_wrap(~lubridate::year(date), scales = "free") +
#     labs(x = expression(R[E]), y = "EWS value", title = do_region) +
#     scale_color_viridis_c(name = "week of year", direction = 1))
# }

No clear signal in those plots. Next I’ll look at the first difference of each EWS (x) between time t and \(t+1\). I calculate the difference in three ways:

  1. Simple difference: \(x(t) - x(t-1)\)
  2. Ratio: \(\frac{x(t)}{x(t-1)}\)
  3. Log ratio: \(\text{log}\left(\frac{x(t)}{x(t-1)}\right)\)
ews_reffs_diffed <- ews_reffs %>%
  group_by(ews, region) %>%
  mutate(
    diff_value = value - lag(value, 1),
    ratio_value = value / lag(value, 1),
    log_ratio_value = log(ratio_value)
  ) %>%
  rename("Reff" = med)
NaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs producedNaNs produced

Correlations between \(R_E\) and first differenced EWS

ews_reffs_diffed %>%
  filter(Reff < 1) %>%
  ggplot(aes(x = Reff, y = round(diff_value,1))) +
  geom_point(aes(color = lubridate::week(date)), alpha = 0.7) +
  facet_wrap(ews~region, scales = "free", ncol = 4) +
  scale_color_viridis_c(name = "Week of year", limits = c(1,52), 
                        breaks = c(1,26,52)) +
  labs(y = "First difference of EWS", x = expression(R[E])) +
  theme(legend.position = "top")

Correlations between \(R_E\) and “first ratioed” EWS

ews_reffs_diffed %>%
  filter(Reff < 1) %>%
  ggplot(aes(x = Reff, y = round(ratio_value,1))) +
  geom_point(aes(color = lubridate::week(date)), alpha = 0.7) +
  facet_wrap(ews~region, scales = "free", ncol = 4) +
  scale_color_viridis_c(name = "Week of year", limits = c(1,52), 
                        breaks = c(1,26,52)) +
  labs(y = "First ratio of EWS", x = expression(R[E])) +
  theme(legend.position = "top")

Correlations between \(R_E\) and log of “first ratioed” EWS

ews_reffs_diffed %>%
  filter(Reff < 1) %>%
  ggplot(aes(x = Reff, y = log_ratio_value)) +
  geom_point(aes(color = lubridate::week(date)), alpha = 0.7) +
  facet_wrap(ews~region, scales = "free", ncol = 4) +
  scale_color_viridis_c(name = "Week of year", limits = c(1,52), 
                        breaks = c(1,26,52)) +
  labs(y = "log(First ratio of EWS)", x = expression(R[E])) +
  theme(legend.position = "top")

Correlations between \(R_E\) and log of “first ratioed” EWS, near-zeros removed

It looks like there might be trends when the EWS are actually “moving”, that is, when \(\text{log}\left(\frac{x(t)}{x(t-1)} \right) \ne 0\) or approximately so. So, I’ll remove those values that are close to zero by assuming all values less than \(|0.05|\) are essentially zero. Now some patterns emerge as expected.

ews_reffs_diffed %>%
  filter(Reff < 1) %>%
  filter(round(log_ratio_value,1) != 0) %>%
  ggplot(aes(x = Reff, y = log_ratio_value)) +
  geom_point(aes(color = lubridate::week(date)), alpha = 0.7) +
  stat_smooth(method = "lm", se = FALSE, linetype = 1, color = "black") +
  facet_wrap(ews~region, scales = "free", ncol = 4) +
  scale_color_viridis_c(name = "Week of year", limits = c(1,52), 
                        breaks = c(1,26,52)) +
  labs(y = "log(First ratio of EWS)", x = expression(R[E])) +
  theme(legend.position = "top")

Look at lag in dynamics

Due to the stochastic seasonal dynamics, there may be a lag between \(R_E\) and EWS that are calculated from raw case data. To look at this, I will plot \(R_E\) and each EWS on the same plot, with peaks in each year highlighted. These plots will come from simulated series where we know that the \(R_E\) and the dynamics are explicitly linked (unlike \(R_E\) estimates from particle filtering with the real data).

all_sim_files <- list.files("../../simulations/")
sim_ids <- grep("bootstrap*", all_sim_files)
sim_files <- all_sim_files[sim_ids]
all_sims <- tibble()  # empty tibble for storage
for(do_file in sim_files){
  tmp_sim <- readRDS(paste0("../../simulations/", do_file)) %>%
    filter(sim == 1) %>%
    unnest() %>%
    mutate(city = str_sub(do_file, 16, -5))
  
  all_sims <- bind_rows(all_sims, tmp_sim)
}
sim_ews <- tibble()  # empty storage df
for(do_city in unique(all_sims$city)){
  tmp_data <- all_sims %>%
    filter(city == do_city)
  
  tmp_cases <- tmp_data$reports
  
  tmp_stats <- spaero::get_stats(
    x = tmp_cases,
    center_trend = "local_constant", 
    center_kernel = "uniform", 
    center_bandwidth = 30, 
    stat_trend = "local_constant", 
    stat_kernel = "uniform", 
    stat_bandwidth = 30, 
    lag = 1, 
    backward_only = TRUE
  )$stats
  
  sim_ews <- sim_ews %>%
    bind_rows(
      tmp_data %>% bind_cols(as_tibble(tmp_stats))
    )
}
sim_ews <- sim_ews %>%
  mutate(critical = ifelse(RE_seas < 1, FALSE, TRUE))
for(do_city in unique(sim_ews$city)){
  tmp <- filter(sim_ews, city == do_city)
  x <- tmp$time
  y <- tmp$RE_seas
  y2 <- log(tmp$variance)
  color <- ifelse(tmp$critical == TRUE, "red", "blue")
  
  par(mar = c(5,5,2,5))
  plot(x, y, type = "n", ylab = expression(italic(R)[E]), xlab = "Date", main = do_city)
  segments(x0 = x, y0 = y, x1 = x+diff(x), y1 = y+diff(y), col = color)
  par(new = T)
  plot(x, y2, type = "l", axes=F, xlab=NA, ylab=NA, lwd = 1.5)
  axis(side = 4)
  mtext(side = 4, line = 3, "Log variance")
}

longer object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object length

longer object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object length

longer object length is not a multiple of shorter object lengthlonger object length is not a multiple of shorter object length

LS0tCnRpdGxlOiAiQ29ycmVsYXRpb25zIGJldHdlZW4gRVdTIGFuZCBlZmZlY3RpdmUgcmVwcm9kdWN0aW9uIG51bWJlciIKb3V0cHV0OiBodG1sX25vdGVib29rCmRhdGU6IDIwMTktMDEtMjIKLS0tCgpUbyBjb21wbGVtZW50IHRoZSBmaXhlZCB3aW5kb3cgYW5hbHlzaXMgZm9jdXNlZCBvbiBlbWVyZ2VuY2UgYW5kIGVsaW1pbmF0aW9uIGV2ZW50cywgd2UgYWxzbyB3YW50IHRvIHNob3cgdGhhdCB0aGUgZHluYW1pY3Mgb2YgRVdTIHJlZmxlY3QgdGhvc2Ugb2YgJFJfRSQuCkkgZXN0aW1hdGVkICRSX0UkIG92ZXIgdGhlIG9ic2VydmF0aW9ucyB2aWEgcGFydGljbGUgZmlsdGVyaW5nIGFuZCBzdGF0ZSBlc3RpbWF0aW9uLgpIZXJlLCBJIHJlYWQgaW4gdGhvc2UgdmFsdWVzIGZvciBlYWNoIGNpdHkgYW5kIGNvcnJlbGF0ZSB0aGVtIHdpdGggdG9wLXBlcmZvcm1pbmcgRVdTLgpUb3AtcGVyZm9ybWluZyBFV1MgYXJlIHRob3NlIHdpdGggaGlnaCBBVUMgaW4gdGhlIGVtZXJnZW5jZSBhbmQgZWxpbWluYXRpb24gc2NlbmFyaW9zLgpOb3RlIHRoYXQgJFJfRSQgd2FzIGNhbGN1bGF0ZWQgYXM6CgokJApSX0UgPSBcZnJhY3tcYmV0YV90fXtcZ2FtbWF9IFxmcmFje1NfdH17Tl90fS4KJCQKCmBgYHtyIGxvYWQtcmV9CmxpYnJhcnkodGlkeXZlcnNlKQoKY2l0aWVzIDwtIGMoIkFnYWRleiIsICJNYXJhZGkiLCAiTmlhbWV5IiwgIlppbmRlciIpICAjIHRoZSBmb3VyIGNpdHkgbmFtZXMKcmVmZnMgPC0gdGliYmxlKCkgICMgZW1wdHkgdGliYmxlIGZvciBzdG9yYWdlCgpmb3IoZG9fY2l0eSBpbiBjaXRpZXMpewogIGZuYW1lIDwtIHBhc3RlMCgiLi4vLi4vcmVzdWx0cy9maWx0ZXJlZC1zdGF0ZXMtIiwgZG9fY2l0eSwgIi5SRFMiKQogIHJhd19maWx0ZXIgPC0gcmVhZFJEUyhmbmFtZSkKICByZWZmX2lkIDwtIHdoaWNoKHJhd19maWx0ZXIkc3RhdGUgPT0gImVmZmVjdGl2ZV9yX3NlYXNvbmFsIikKICB0bXBfcmVmZiA8LSByYXdfZmlsdGVyJGRhdGFbW3JlZmZfaWRdXSAlPiUKICAgIGRwbHlyOjpzZWxlY3QodGltZSwgbWVkLCB3ZWVrLCBkYXRlKSAlPiUKICAgIGRwbHlyOjptdXRhdGUoY2l0eSA9IGRvX2NpdHkpCiAgCiAgcmVmZnMgPC0gYmluZF9yb3dzKHJlZmZzLCB0bXBfcmVmZikKfQoKZ2dwbG90KHJlZmZzLCBhZXMoeCA9IGRhdGUsIHkgPSBtZWQsIGNvbG9yID0gY2l0eSkpICsKICBnZW9tX2hsaW5lKGFlcyh5aW50ZXJjZXB0ID0gMSksIGxpbmV0eXBlID0gMiwgY29sb3IgPSAiZ3JleTM1IikgKwogIGdlb21fbGluZSgpICsKICBsYWJzKHggPSAiZGF0ZSIsIHkgPSBleHByZXNzaW9uKGl0YWxpYyhSKVtFXSksIAogICAgICAgdGl0bGUgPSAiRWZmZWN0aXZlIHJlcHJvZHVjdGlvbiBudW1iZXIgb3ZlciB0aW1lIikgKwogIGdndGhlbWVzOjpzY2FsZV9jb2xvcl9jb2xvcmJsaW5kKCkKYGBgCgpOb3cgSSBjYWxjdWxhdGUgdGhlIEVXUyBhdCBlYWNoIG9ic2VydmF0aW9uIHRpbWUgdXNpbmcgdGhlIGBzcGFlcm86OmdldF9zdGF0cygpYCBmdW5jdGlvbi4KSSB1c2UgYSBiYW5kd2lkdGggb2YgMzAgd2Vla3MgYW5kIHNldCBgYmFja3dhcmRfb25seSA9IFRSVUVgIHNvIHRoYXQgb25seSBwYXN0IGRhdGEgYXJlIHVzZWQgd2hlbiBjYWxjdWxhdGluZyBFV1MuCgpgYGB7ciBld3MtY2FsY30KbGlicmFyeShzcGFlcm8pCgpmbmFtZSA8LSAiLi4vLi4vZGF0YS9jbGVhbi1kYXRhL3dlZWtseS1tZWFzbGVzLWluY2lkZW5jZS1uaWdlci1jaXRpZXMtY2xlYW4uUkRTIgptZWFzbGVzX2RhdGEgPC0gcmVhZFJEUyhmbmFtZSkgJT4lCiAgZmlsdGVyKHllYXIgPiAxOTk0KSAgIyBkcm9wIGZpcnN0IE5BIHllYXIsIG9ubHkgdXNlZCBmb3IgbW9kZWxpbmcKCmFsbF9zdGF0cyA8LSB0aWJibGUoKSAgIyBlbXB0eSB0aWJibGUgZm9yIHN0b3JhZ2UKCmZvcihkb19yZWdpb24gaW4gdW5pcXVlKG1lYXNsZXNfZGF0YSRyZWdpb24pKXsKICAKICBjYXNlcyA8LSBtZWFzbGVzX2RhdGEgJT4lCiAgICBmaWx0ZXIocmVnaW9uID09IGRvX3JlZ2lvbikgJT4lCiAgICBwdWxsKGNhc2VzKQogIAogIGNpdHlfc3RhdHMgPC0gc3BhZXJvOjpnZXRfc3RhdHMoCiAgICB4ID0gY2FzZXMsCiAgICBjZW50ZXJfdHJlbmQgPSAibG9jYWxfY29uc3RhbnQiLCAKICAgIGNlbnRlcl9rZXJuZWwgPSAidW5pZm9ybSIsIAogICAgY2VudGVyX2JhbmR3aWR0aCA9IDMwLCAKICAgIHN0YXRfdHJlbmQgPSAibG9jYWxfY29uc3RhbnQiLCAKICAgIHN0YXRfa2VybmVsID0gInVuaWZvcm0iLCAKICAgIHN0YXRfYmFuZHdpZHRoID0gMzAsIAogICAgbGFnID0gMSwgCiAgICBiYWNrd2FyZF9vbmx5ID0gVFJVRQogICkkc3RhdHMKICAKICBjaXR5X3N0YXRzX3RiIDwtIGFzX3RpYmJsZShjaXR5X3N0YXRzKSAlPiUKICAgIG11dGF0ZSgKICAgICAgdGltZV9pdGVyID0gMTpuKCksCiAgICAgIGRhdGUgPSB1bmlxdWUobWVhc2xlc19kYXRhJGRhdGUpCiAgICApICU+JQogICAgZ2F0aGVyKGtleSA9IGV3cywgdmFsdWUgPSB2YWx1ZSwgLXRpbWVfaXRlciwgLWRhdGUpICU+JQogICAgbXV0YXRlKHJlZ2lvbiA9IGRvX3JlZ2lvbikKICAKICBhbGxfc3RhdHMgPC0gYmluZF9yb3dzKGFsbF9zdGF0cywgY2l0eV9zdGF0c190YikKfQoKIyBiZXN0X2V3cyA8LSBjKCJhdXRvY29ycmVsYXRpb24iLCAiYXV0b2NvdmFyaWFuY2UiLCAibWVhbiIsICJ2YXJpYW5jZSIpCiMgYWxsX3N0YXRzIDwtIGFsbF9zdGF0cyAlPiUKIyAgIGRwbHlyOjpmaWx0ZXIoZXdzICVpbiUgYmVzdF9ld3MpCmBgYAoKTm93IEkgbG9vayBhdCBzY2F0dGVycGxvdHMgb2YgJFJfRSQgdmVyc3VzIEVXUy4KCmBgYHtyIHNjYXR0ZXJzfQpyZWZmcyA8LSByZWZmcyAlPiUKICBkcGx5cjo6bXV0YXRlKHJlZ2lvbiA9IHBhc3RlMChjaXR5LCAiIChDaXR5KSIpKSAlPiUKICBkcGx5cjo6c2VsZWN0KC1jaXR5KQoKZXdzX3JlZmZzIDwtIGxlZnRfam9pbihhbGxfc3RhdHMsIHJlZmZzKQpld3NfcmVmZnNfc3ViY3JpdCA8LSBld3NfcmVmZnMgJT4lCiAgZmlsdGVyKG1lZCA8IDEpCgpnZ3Bsb3QoZXdzX3JlZmZzX3N1YmNyaXQsIGFlcyh4ID0gbWVkLCB5ID0gdmFsdWUsIGNvbG9yID0gcmVnaW9uKSkgKwogIGdlb21fcG9pbnQoYWxwaGEgPSAwLjIpICsKICBmYWNldF93cmFwKH5ld3MsIHNjYWxlcyA9ICJmcmVlIikgKwogIGdndGhlbWVzOjpzY2FsZV9jb2xvcl9jb2xvcmJsaW5kKCkKCiMgZm9yKGRvX3JlZ2lvbiBpbiB1bmlxdWUoZXdzX3JlZmZzJHJlZ2lvbikpewojICAgcHJpbnQoZ2dwbG90KGZpbHRlcihld3NfcmVmZnMsIHJlZ2lvbiA9PSBkb19yZWdpb24gJiBld3MgPT0gInZhcmlhbmNlIiksIAojICAgICAgICBhZXMoeCA9IG1lZCwgeSA9IHZhbHVlLCBjb2xvciA9IGx1YnJpZGF0ZTo6d2VlayhkYXRlKSkpICsKIyAgICAgZ2VvbV9wb2ludCgpICsKIyAgICAgZmFjZXRfd3JhcCh+bHVicmlkYXRlOjp5ZWFyKGRhdGUpLCBzY2FsZXMgPSAiZnJlZSIpICsKIyAgICAgbGFicyh4ID0gZXhwcmVzc2lvbihSW0VdKSwgeSA9ICJFV1MgdmFsdWUiLCB0aXRsZSA9IGRvX3JlZ2lvbikgKwojICAgICBzY2FsZV9jb2xvcl92aXJpZGlzX2MobmFtZSA9ICJ3ZWVrIG9mIHllYXIiLCBkaXJlY3Rpb24gPSAxKSkKIyB9CmBgYAoKTm8gY2xlYXIgc2lnbmFsIGluIHRob3NlIHBsb3RzLgpOZXh0IEknbGwgbG9vayBhdCB0aGUgZmlyc3QgZGlmZmVyZW5jZSBvZiBlYWNoIEVXUyAoKngqKSBiZXR3ZWVuIHRpbWUgKnQqIGFuZCAkdCsxJC4KSSBjYWxjdWxhdGUgdGhlIGRpZmZlcmVuY2UgaW4gdGhyZWUgd2F5czoKCjEuIFNpbXBsZSBkaWZmZXJlbmNlOiAkeCh0KSAtIHgodC0xKSQKMi4gUmF0aW86ICRcZnJhY3t4KHQpfXt4KHQtMSl9JAozLiBMb2cgcmF0aW86ICRcdGV4dHtsb2d9XGxlZnQoXGZyYWN7eCh0KX17eCh0LTEpfVxyaWdodCkkCgpgYGB7ciBjYWxjLWV3cy1kaWZmc30KZXdzX3JlZmZzX2RpZmZlZCA8LSBld3NfcmVmZnMgJT4lCiAgZ3JvdXBfYnkoZXdzLCByZWdpb24pICU+JQogIG11dGF0ZSgKICAgIGRpZmZfdmFsdWUgPSB2YWx1ZSAtIGxhZyh2YWx1ZSwgMSksCiAgICByYXRpb192YWx1ZSA9IHZhbHVlIC8gbGFnKHZhbHVlLCAxKSwKICAgIGxvZ19yYXRpb192YWx1ZSA9IGxvZyhyYXRpb192YWx1ZSkKICApICU+JQogIHJlbmFtZSgiUmVmZiIgPSBtZWQpCmBgYAoKIyMjIENvcnJlbGF0aW9ucyBiZXR3ZWVuICRSX0UkIGFuZCBmaXJzdCBkaWZmZXJlbmNlZCBFV1MKCmBgYHtyIGZpcnN0LWRpZmYtc2NhdHRlcnMsIGZpZy5oZWlnaHQ9MjAsIGZpZy53aWR0aD04fQpld3NfcmVmZnNfZGlmZmVkICU+JQogIGZpbHRlcihSZWZmIDwgMSkgJT4lCiAgZ2dwbG90KGFlcyh4ID0gUmVmZiwgeSA9IHJvdW5kKGRpZmZfdmFsdWUsMSkpKSArCiAgZ2VvbV9wb2ludChhZXMoY29sb3IgPSBsdWJyaWRhdGU6OndlZWsoZGF0ZSkpLCBhbHBoYSA9IDAuNykgKwogIGZhY2V0X3dyYXAoZXdzfnJlZ2lvbiwgc2NhbGVzID0gImZyZWUiLCBuY29sID0gNCkgKwogIHNjYWxlX2NvbG9yX3ZpcmlkaXNfYyhuYW1lID0gIldlZWsgb2YgeWVhciIsIGxpbWl0cyA9IGMoMSw1MiksIAogICAgICAgICAgICAgICAgICAgICAgICBicmVha3MgPSBjKDEsMjYsNTIpKSArCiAgbGFicyh5ID0gIkZpcnN0IGRpZmZlcmVuY2Ugb2YgRVdTIiwgeCA9IGV4cHJlc3Npb24oUltFXSkpICsKICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAidG9wIikKYGBgCgojIyMgQ29ycmVsYXRpb25zIGJldHdlZW4gJFJfRSQgYW5kICJmaXJzdCByYXRpb2VkIiBFV1MKCmBgYHtyIGZpZy5oZWlnaHQ9MjAsIGZpZy53aWR0aD04fQpld3NfcmVmZnNfZGlmZmVkICU+JQogIGZpbHRlcihSZWZmIDwgMSkgJT4lCiAgZ2dwbG90KGFlcyh4ID0gUmVmZiwgeSA9IHJvdW5kKHJhdGlvX3ZhbHVlLDEpKSkgKwogIGdlb21fcG9pbnQoYWVzKGNvbG9yID0gbHVicmlkYXRlOjp3ZWVrKGRhdGUpKSwgYWxwaGEgPSAwLjcpICsKICBmYWNldF93cmFwKGV3c35yZWdpb24sIHNjYWxlcyA9ICJmcmVlIiwgbmNvbCA9IDQpICsKICBzY2FsZV9jb2xvcl92aXJpZGlzX2MobmFtZSA9ICJXZWVrIG9mIHllYXIiLCBsaW1pdHMgPSBjKDEsNTIpLCAKICAgICAgICAgICAgICAgICAgICAgICAgYnJlYWtzID0gYygxLDI2LDUyKSkgKwogIGxhYnMoeSA9ICJGaXJzdCByYXRpbyBvZiBFV1MiLCB4ID0gZXhwcmVzc2lvbihSW0VdKSkgKwogIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICJ0b3AiKQpgYGAKCiMjIyBDb3JyZWxhdGlvbnMgYmV0d2VlbiAkUl9FJCBhbmQgbG9nIG9mICJmaXJzdCByYXRpb2VkIiBFV1MKCmBgYHtyIGZpZy5oZWlnaHQ9MjAsIGZpZy53aWR0aD04fQpld3NfcmVmZnNfZGlmZmVkICU+JQogIGZpbHRlcihSZWZmIDwgMSkgJT4lCiAgZ2dwbG90KGFlcyh4ID0gUmVmZiwgeSA9IGxvZ19yYXRpb192YWx1ZSkpICsKICBnZW9tX3BvaW50KGFlcyhjb2xvciA9IGx1YnJpZGF0ZTo6d2VlayhkYXRlKSksIGFscGhhID0gMC43KSArCiAgZmFjZXRfd3JhcChld3N+cmVnaW9uLCBzY2FsZXMgPSAiZnJlZSIsIG5jb2wgPSA0KSArCiAgc2NhbGVfY29sb3JfdmlyaWRpc19jKG5hbWUgPSAiV2VlayBvZiB5ZWFyIiwgbGltaXRzID0gYygxLDUyKSwgCiAgICAgICAgICAgICAgICAgICAgICAgIGJyZWFrcyA9IGMoMSwyNiw1MikpICsKICBsYWJzKHkgPSAibG9nKEZpcnN0IHJhdGlvIG9mIEVXUykiLCB4ID0gZXhwcmVzc2lvbihSW0VdKSkgKwogIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICJ0b3AiKQpgYGAKCiMjIyBDb3JyZWxhdGlvbnMgYmV0d2VlbiAkUl9FJCBhbmQgbG9nIG9mICJmaXJzdCByYXRpb2VkIiBFV1MsIG5lYXItemVyb3MgcmVtb3ZlZAoKSXQgbG9va3MgbGlrZSB0aGVyZSBtaWdodCBiZSB0cmVuZHMgd2hlbiB0aGUgRVdTIGFyZSBhY3R1YWxseSAibW92aW5nIiwgdGhhdCBpcywgd2hlbiAkXHRleHR7bG9nfVxsZWZ0KFxmcmFje3godCl9e3godC0xKX0gXHJpZ2h0KSBcbmUgMCQgb3IgYXBwcm94aW1hdGVseSBzby4KU28sIEknbGwgcmVtb3ZlIHRob3NlIHZhbHVlcyB0aGF0IGFyZSBjbG9zZSB0byB6ZXJvIGJ5IGFzc3VtaW5nIGFsbCB2YWx1ZXMgbGVzcyB0aGFuICR8MC4wNXwkIGFyZSBlc3NlbnRpYWxseSB6ZXJvLgpOb3cgc29tZSBwYXR0ZXJucyBlbWVyZ2UgYXMgZXhwZWN0ZWQuCgpgYGB7ciBmaWcuaGVpZ2h0PTIwLCBmaWcud2lkdGg9OH0KZXdzX3JlZmZzX2RpZmZlZCAlPiUKICBmaWx0ZXIoUmVmZiA8IDEpICU+JQogIGZpbHRlcihyb3VuZChsb2dfcmF0aW9fdmFsdWUsMSkgIT0gMCkgJT4lCiAgZ2dwbG90KGFlcyh4ID0gUmVmZiwgeSA9IGxvZ19yYXRpb192YWx1ZSkpICsKICBnZW9tX3BvaW50KGFlcyhjb2xvciA9IGx1YnJpZGF0ZTo6d2VlayhkYXRlKSksIGFscGhhID0gMC43KSArCiAgc3RhdF9zbW9vdGgobWV0aG9kID0gImxtIiwgc2UgPSBGQUxTRSwgbGluZXR5cGUgPSAxLCBjb2xvciA9ICJibGFjayIpICsKICBmYWNldF93cmFwKGV3c35yZWdpb24sIHNjYWxlcyA9ICJmcmVlIiwgbmNvbCA9IDQpICsKICBzY2FsZV9jb2xvcl92aXJpZGlzX2MobmFtZSA9ICJXZWVrIG9mIHllYXIiLCBsaW1pdHMgPSBjKDEsNTIpLCAKICAgICAgICAgICAgICAgICAgICAgICAgYnJlYWtzID0gYygxLDI2LDUyKSkgKwogIGxhYnMoeSA9ICJsb2coRmlyc3QgcmF0aW8gb2YgRVdTKSIsIHggPSBleHByZXNzaW9uKFJbRV0pKSArCiAgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gInRvcCIpCmBgYAoKIyBMb29rIGF0IGxhZyBpbiBkeW5hbWljcwpEdWUgdG8gdGhlIHN0b2NoYXN0aWMgc2Vhc29uYWwgZHluYW1pY3MsIHRoZXJlIG1heSBiZSBhIGxhZyBiZXR3ZWVuICRSX0UkIGFuZCBFV1MgdGhhdCBhcmUgY2FsY3VsYXRlZCBmcm9tIHJhdyBjYXNlIGRhdGEuClRvIGxvb2sgYXQgdGhpcywgSSB3aWxsIHBsb3QgJFJfRSQgYW5kIGVhY2ggRVdTIG9uIHRoZSBzYW1lIHBsb3QsIHdpdGggcGVha3MgaW4gZWFjaCB5ZWFyIGhpZ2hsaWdodGVkLgpUaGVzZSBwbG90cyB3aWxsIGNvbWUgZnJvbSBzaW11bGF0ZWQgc2VyaWVzIHdoZXJlIHdlIGtub3cgdGhhdCB0aGUgJFJfRSQgYW5kIHRoZSBkeW5hbWljcyBhcmUgZXhwbGljaXRseSBsaW5rZWQgKHVubGlrZSAkUl9FJCBlc3RpbWF0ZXMgZnJvbSBwYXJ0aWNsZSBmaWx0ZXJpbmcgd2l0aCB0aGUgcmVhbCBkYXRhKS4KCmBgYHtyIGxvYWQtc2ltc30KYWxsX3NpbV9maWxlcyA8LSBsaXN0LmZpbGVzKCIuLi8uLi9zaW11bGF0aW9ucy8iKQpzaW1faWRzIDwtIGdyZXAoImJvb3RzdHJhcCoiLCBhbGxfc2ltX2ZpbGVzKQpzaW1fZmlsZXMgPC0gYWxsX3NpbV9maWxlc1tzaW1faWRzXQoKYWxsX3NpbXMgPC0gdGliYmxlKCkgICMgZW1wdHkgdGliYmxlIGZvciBzdG9yYWdlCmZvcihkb19maWxlIGluIHNpbV9maWxlcyl7CiAgdG1wX3NpbSA8LSByZWFkUkRTKHBhc3RlMCgiLi4vLi4vc2ltdWxhdGlvbnMvIiwgZG9fZmlsZSkpICU+JQogICAgZmlsdGVyKHNpbSA9PSAxKSAlPiUKICAgIHVubmVzdCgpICU+JQogICAgbXV0YXRlKGNpdHkgPSBzdHJfc3ViKGRvX2ZpbGUsIDE2LCAtNSkpCiAgCiAgYWxsX3NpbXMgPC0gYmluZF9yb3dzKGFsbF9zaW1zLCB0bXBfc2ltKQp9CmBgYAoKYGBge3IgY2FsYy1zaW0tZXdzfQpzaW1fZXdzIDwtIHRpYmJsZSgpICAjIGVtcHR5IHN0b3JhZ2UgZGYKZm9yKGRvX2NpdHkgaW4gdW5pcXVlKGFsbF9zaW1zJGNpdHkpKXsKICB0bXBfZGF0YSA8LSBhbGxfc2ltcyAlPiUKICAgIGZpbHRlcihjaXR5ID09IGRvX2NpdHkpCiAgCiAgdG1wX2Nhc2VzIDwtIHRtcF9kYXRhJHJlcG9ydHMKICAKICB0bXBfc3RhdHMgPC0gc3BhZXJvOjpnZXRfc3RhdHMoCiAgICB4ID0gdG1wX2Nhc2VzLAogICAgY2VudGVyX3RyZW5kID0gImxvY2FsX2NvbnN0YW50IiwgCiAgICBjZW50ZXJfa2VybmVsID0gInVuaWZvcm0iLCAKICAgIGNlbnRlcl9iYW5kd2lkdGggPSAzMCwgCiAgICBzdGF0X3RyZW5kID0gImxvY2FsX2NvbnN0YW50IiwgCiAgICBzdGF0X2tlcm5lbCA9ICJ1bmlmb3JtIiwgCiAgICBzdGF0X2JhbmR3aWR0aCA9IDMwLCAKICAgIGxhZyA9IDEsIAogICAgYmFja3dhcmRfb25seSA9IFRSVUUKICApJHN0YXRzCiAgCiAgc2ltX2V3cyA8LSBzaW1fZXdzICU+JQogICAgYmluZF9yb3dzKAogICAgICB0bXBfZGF0YSAlPiUgYmluZF9jb2xzKGFzX3RpYmJsZSh0bXBfc3RhdHMpKQogICAgKQp9CmBgYAoKYGBge3IgcGxvdC1keW5hbWljcywgZmlnLmhlaWdodD00LCBmaWcud2lkdGg9MTB9CnNpbV9ld3MgPC0gc2ltX2V3cyAlPiUKICBtdXRhdGUoY3JpdGljYWwgPSBpZmVsc2UoUkVfc2VhcyA8IDEsIEZBTFNFLCBUUlVFKSkKCmZvcihkb19jaXR5IGluIHVuaXF1ZShzaW1fZXdzJGNpdHkpKXsKICB0bXAgPC0gZmlsdGVyKHNpbV9ld3MsIGNpdHkgPT0gZG9fY2l0eSkKICB4IDwtIHRtcCR0aW1lCiAgeSA8LSB0bXAkUkVfc2VhcwogIHkyIDwtIGxvZyh0bXAkdmFyaWFuY2UpCiAgY29sb3IgPC0gaWZlbHNlKHRtcCRjcml0aWNhbCA9PSBUUlVFLCAicmVkIiwgImJsdWUiKQogIAogIHBhcihtYXIgPSBjKDUsNSwyLDUpKQogIHBsb3QoeCwgeSwgdHlwZSA9ICJuIiwgeWxhYiA9IGV4cHJlc3Npb24oaXRhbGljKFIpW0VdKSwgeGxhYiA9ICJEYXRlIiwgbWFpbiA9IGRvX2NpdHkpCiAgc2VnbWVudHMoeDAgPSB4LCB5MCA9IHksIHgxID0geCtkaWZmKHgpLCB5MSA9IHkrZGlmZih5KSwgY29sID0gY29sb3IpCiAgcGFyKG5ldyA9IFQpCiAgcGxvdCh4LCB5MiwgdHlwZSA9ICJsIiwgYXhlcz1GLCB4bGFiPU5BLCB5bGFiPU5BLCBsd2QgPSAxLjUpCiAgYXhpcyhzaWRlID0gNCkKICBtdGV4dChzaWRlID0gNCwgbGluZSA9IDMsICJMb2cgdmFyaWFuY2UiKQp9CgoKYGBg